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ABSTRACT 


Established,  steady-state  flow  beyond  the  entry  region  in  a  duct  with 
periodic  wall  geometry  is  hypothesized  as  spatially  periodic.  Spatially 
periodic  flow  is  considered  to  be  attained  when  the  flow  and  heat  transfer 
in  any  cycle  of  periodic  wall  geometry  is  identical  to  the  cycle  before 
and  after  it.  Thus,  a  two-dimensional  duct  with  periodic  wall  geometry 
may  have  its  established  flow  far  from  the  entrance  predicted  by  a  numeri¬ 
cal  solution  over  a  cycle  of  the  geometry. 

A  corrugated  duct  and  a  parallel  plate  duct  with  repeated  rib  rough¬ 
ness  were  selected  for  numerical  modeling.  These  duct  configurations  are 
of  interest  in  heat  exchanger  designs  and  other  components  in  propulsion 
systems. 

A  numerical  method,  utilizing  the  general  finite-difference  program 
TEACH  from  Imperial  College,  was  successfully  developed  to  solve  for  fully 
established,  laminar  flow  in  a  corrugated  duct.  This  method  was  sub¬ 
sequently  applied  to  predict  laminar  flow  and  heat  transfer  downstream 
in  a  parallel  plate  duct  with  repeated  ribs  as  surface  roughness  elements. 

Results  for  the  corrugated  duct  reveal  the  presence  of  large  recir¬ 
culation  zones  at  Reynolds  numbers  greater  than  about  50.  The  average 
friction  factor  varies  inversely  with  the  Reynolds  number  at  low  flow 


rates.  It  then  passes  through  a  minimum,  as  the  Reynolds  number  is  varied. 


and  increases,  apparently  approaching  an  approximately  constant  value  coj 


responding  to  the  development  of  the  recirculation  zones. 
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Results  for  the  rough  duct  are  preliminary,  but  indicate  reasonable 
results  may  be  obtained.  Both  the  flow  field  and  the  temperature  dis¬ 
tribution  in  the  solid  plate  are  determined,  thereby  providing  means  to 
predict  the  heat  transfer  behavior  when  the  thermal  resistance  in  the 
fluid  and  the  thermal  resistance  in  the  wall  are  of  the  same  order  of 


magnitude. 
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NOMENCLATURE 


Definition _ 

Distance  between  comer  regions,  a 

Hydraulic  diameter,  a 

Friction  factor 

Corrugated  channel  gap,  a 

Maximum  width  of  corrugated  channel,  a 

Rib  height ,  m 

Thermal  conductivity,  w/n-°C 
Contour  co-ordinate,  a 

Mean  distance  between  corresponding  planes,  a 

Mass  Flow  Rate,  kg/s 
Number  of  nodes  cross-stream 

Nusselt  Number 
2 

Pressure,  N/m 

Average  pressure  drop  between  corresponding  planes 
Peclet  number  (Pe  *  RePr) 

Prandtl  number 

? 

Average  heat  flux  into  cell,  W/nT-s 

2 

Imposed  heat  flux  on  plate,  W/m  -s 

2 

Local  heat  flux  at  fluid-plate  interface,  W/m  -s 
Reynolds  number 


Rib  separation 


Symbol 


S'" 

x 

S'" 

y 


T 


u 

V 

w 

X 
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NOMENCLATURE — Continued 

Definition _ 

3 

Momentum  source  in  x-momentum  equation,  N/m 

3 

Momentum  source  in  y-momentum  equation,  N/m 

3 

Energy  source  in  energy  equation,  W/m 

Temperature,  °C 

Velocity  in  x-direction,  m/s 

Velocity  in  y-direction,  m/s 

Rib  width ,  m 

Cartesian  co-ordinate  direction 

Cartesian  co-ordinate  direction 

Distance  from  centerline  to  plate  wall,  m 

Thickness  of  plate,  m 
3 

Density,  kg/m 
Viscosity,  kg/m-s 

Tolerance  of  core  computer  program 

Overall  computer  program  tolerance 

2 

Kinematic  viscosity,  m  /s 
Resistance  coefficient 


Subscripts 

av  Average  value 

b  Boundary  value 

B  Bulk  value 

f  Fluid  value 


i 


Nodal  index 


NOMENCLATURE— Continued 

Subscripts 

Symbol _ Definition _ 

Internal  plane  value 
Nodal  index 
Plate  value 
Wall  or  local  value 

Superscripts 

Denotes  plane  immediately  upstream  of  corresponding  plane 
Denotes  plane  immediately  downstream  of  corresponding  plane 
Correction  value 
Non-diraensionalized  value 


Estimated  value 


CHAPTER  1 


INTRODUCTION 


When  studying  the  flow  in  ducts,  there  arise  some  problems  for 
which  a  numerical  finite-difference  method  is  the  only  means  of  obtaining 
a  solution,  albeit  approximate.  One  class  of  such  problems,  important 
for  the  conveccive  heat  transfer  in  advanced  concepts  for  shipbound  pro¬ 
pulsion,  concerns  flow  downstream  of  the  entry  region  in  two-dimensional 
ducts  with  repeating  geomecry. 

The  repetition  of  a  geometrical  characteristic  is  outlined  in 
Fig.  1-1.  The  two-dimensional  duct  is  formed  by  a  pair  of  parallel  plates 
with  square  ribs  fastened  to  the  surfaces  perpendicular  to  the  flow.  A 
rib  on  one  surface  opposes  a  counterpart  on  the  ocher  surface. 


k 


Fig.  1-1.  Duct  with  Repeating  Geometry 

This  duct  can  be  divided  into  geometric  units  which  are  identical 
to  one  another.  The  vertical  dashed  lines  separate  the  duct  segment  in 
Fig.  1-1  into  two  identical  geometric  units  which  will  be  referred  to  as 
'cells'  in  the  remainder  of  this  report. 


Established  flow  is  considered  attained  when  the  flow  pattern  in 
any  cell  is  independent  of  the  entry  region  effects  and  identical  to  the 
flow  in  cells  upstream  and  downstream  of  itself.  It  is  hypothesized  the 
repeating  geometry  produces  a  repeated  flow  pattern,  i.e.,  a  spatially 
periodic  flow. 

This  report  deals  with  the  development  and  application  of  a  nu¬ 
merical  technique  incorporating  a  general  computer  program  for  two- 
dimensional,  recirculating  flows  to  produce  a  solution  for  the  flow  in 
such  cells.  Since  the  cell  flow  pattern  may  be  assumed  representative 
of  the  flow  in  the  duct  beyond  the  entrance  region,  ducts  whose  L/D^ 
ratio  is  large  may  have  their  dominant  pressure  drop  and  heat  transfer 
characteristics  predicted  by  a  cell  solution.  Solving  for  the  flow  in 
a  cell  rather  than  a  more  extensive  solution  involving  the  entrance  re¬ 
gion  then  becomes  an  attractive  means  for  obtaining  useful  engineering 
information  at  reduced  computational  costs. 

The  numerical  technique  developed  in  this  report  will  be  applied 
to  two  problems.  The  first  problem  is  concerned  with  flows  in  corru¬ 
gated  channels  and  fractures;  the  second  problem  is  concerned  with  flows 
between  ribbed  parallel  plates.  These  problems  will  be  titled  by  the 
names  of  the  computer  programs  developed  to  solve  them:  CRACFLO  and 
RUFF  respectively. 


1.1  Problem  Statement  for  C RAC FLO 

For  compact  heat  exchangers  using  plate-fin  geometries,  one  means 
of  augmenting  heat  transfer  is  use  of  a  corrugated  surface  for  the  plate 
[Kays  and  London,  1964].  This  geometry  is  sketched  in  Figure  1-2. 

Studies  of  heat  exchangers  with  corrugated  wall  channels  are  few 
in  number.  General  friction  factors  and  Nusselt  numbers  for  a  single 
corrugated  wall  channel  geometry  are  presented  by  Belabarodov  and  Volgin 
[1971].  The  friction  factors  were  determined  for  a  Reynolds  number 
range  of  5000  to  45,000.  The  Nusselt  numbers  were  found  for  4000  <  Re 
45,000.  A  definition  of  the  Reynolds  number  was  not  provided  in  the 
paper. 

Two  other  papers  [Goldstein  and  Sparrow,  1975  and  1977]  investi¬ 
gate  the  mass  and  heat  transfer  in  corrugated  wall  channels.  The  invest¬ 
igations  limit  study  to  the  first  two  cycles  of  the  corrugated  channel 
beyond  the  entrance  and  do  not  provide  pressure  drop  or  friction  factor 
data.  The  range  of  Reynolds  numbers  in  these  studies  was  from  150  to 
8550  where 

2  H  U 

Re  -  — 5  (1-1) 

The  symbol  H  is  the  maximum  channel  width  as  shown  in  Figure  1-2  (c). 

Only  one  geometry  was  investigated. 

Data  for  three  corrugated  wall  geometries  are  presented  by  Kays 
and  London,  [1964],  Average  friction  factors  and  Stanton- [Prandtl]  2/3 


)  Plate-fin  heat  exchanqer  qeometry 
with  corrugated  surface 


b)  Conceptualized  fracture 


c)  Corrugated  wall  channel 


L-2.  Conceptualized  Fracture  and  Corrugated 
Wall  Channel  Schematics 


number  products  are  presented  for  the  Reynolds  number  range  of  300  to 
10,000.  The  Reynolds  number  is  based  on  a  hydraulic  diameter  defined  by: 


4  (minimum  flow  area)  (streamwise  length) 
(transfer  surface  area) 


(1-2) 


This  brief  review  represents  the  sum  of  the  literature  on  corru¬ 
gated  wall  channels  located  by  the  author.  Very  few  geometries  have  been 
studied  and  the  data  are  incomplete  as  average  friction  factors  have  not 
even  been  presented  for  laminar  flow  in  all  cases. 

A  related  problem  is  the  flow  in  fractures.  There  arises,  on 
occasion,  a  need  or  desire  for  average  friction  factor  versus  Reynolds 
number  or  average  friction  factor  versus  geometry  relationships  for  frac¬ 
tures  and  cracks.  One  such  case  is  for  cracks  which  may  occur  in  pres¬ 
surized  components  of  power  systems  [Button  et  al.,  1978]. 

A  crack  in  a  pressurized  component  is  a  warning  of  impending 
failure.  If  gases  escaping  through  this  crack  may  be  detected,  an  early 
warning  of  the  component's  failure  is  available.  The  construction  of  a 
leak  detector  to  serve  as  a  warning  device  requires  an  estimate  to  be 
made  of  the  leak  race  of  gases  through  a  crack.  The  leak  rate  is  depen¬ 
dent  on,  among  other  things,  the  average  friction  factor. 

Flow  in  rock  fractures  is  another  case  where  the  prediction  of 
Che  average  friction  factor  is  desirable.  The  Dry-Rock  geothermal  pro¬ 
ject  is  a  means  for  extracting  energy  by  circulating  water  in  a  fracture 
system  in  hot,  dry  rock  deep  below  the  surface  of  the  earth  [McFarland, 
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1975].  The  overall  flow  impedance  of  the  fracture  system  is  important 
since  it  determines  the  pumping  power  required  to  flow  water  through 
the  fracture  system.  The  local  fracture  impedance  has  importance  when 
investigating  the  circulation  patterns  of  the  flow. 

Two  prediction  methods  available  [McFarland,  1975;  McFarland  and 
Murphy,  1976]  use  the  D'Arcy  relation  for  hydraulic  diffusion  to  simplify 
the  governing  equations.  The  D'Arcy  relation  is  based  on  a  friction 
factor-Reynolds  number  relationship  similar  to  that  for  laminar  flow 
between  parallel  plates.  Nemat-Nasser  and  Ohtsubo  [1978]  assumed  a  lam¬ 
inar  parabolic  profile  to  occur  across  the  fracture  gap,  neglecting  the 
rough  character  of  the  surface. 

Therefore,  current  predictions  of  flow  patterns  in  cracks  and, 
consequently,  the  extraction  of  thermal  energy  from  them  are  based  on  a 
simplified  model  of  the  local  fracture  impedance.  It  is  desirable  to 
determine  the  local  fracture  impedance  more  realistically  to  improve  the 
accuracy  of  predictions. 

The  local  flow  in  a  fracture  may  be  assumed  to  exhibit  two- 
dimensional  behavior  in  regions  suitably  far  removed  from  the  inlet  or 
outlet.  Fractures  do  not  exhibit  repeating  geometry,  but  have  the  poten¬ 
tial  for  their  average  characteristics  to  be  approximated  by  a  repeating 
geometry.  A  conceptualized  fracture  and  corrugated  wall  channel  are 
shown  in  Figure  1-2.  An  approximation  of  the  fracture  may  be  achieved 
through  careful  selection  of  a  corrugated  wall  channel  geometry.  The 
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average  characteristics  of  a  crack  or  fracture  may  then  be  estimated  by 
performing  a  numerical  solution  of  the  flow  in  a  cell  of  the  corrugated 
wall  channel. 

Therefore,  it  is  desirable  to  obtain  an  understanding  of  the 
flow  of  fluids  within  a  corrugated  wall  channel.  Furthermore,  since  flow 
in  compact  heat  exchangers  and  fractures  is  often  in  the  laminar  flow 
regime,  a  solution  for  laminar  flow  is  sought.  Program  CRACFLO  has  been 
developed  to  provide  a  solution  for  established,  laminar  flow  in  cor¬ 
rugated  wall  channels  and,  in  particular,  to  provide  average  friction 
factors  based  on  calculated  pressure  drops. 

Program  CRACFLO  utilizes  two  geometric  units  or  cells.  These 
cells  are  shown  in  Figure  1-2  (b)  delimited  by  the  dashed  lines.  The 
two  cell  section  is  shown  rotated  45  degrees  clockwise  in  Figure  1-3. 
Characteristic  dimensions  given  in  this  figure  are  g,  the  channel  gap 
width,  and  a,  the  distance  between  the  corner  regions.  The  geometric 
restrictions  for  this  problem  are:  1)  the  corrugations  must  include 
90  degrees,  i.e.,  square  corners  and  2)  the  distance  between  the  comer 
regions,  a,  cannot  be  less  than  zero. 

To  enable  modeling  of  many  possible  fracture  geometries,  the 
length  a  is  variable  while  g  is  held  constant.  The  geometry  is  then 
characterized  by  the  ratio  a/g.  The  flow  is  characterized  by  a  Rey¬ 
nolds  number 


Re 


p  uB  g 


u 

where  UB  represents  the  bulk  velocity  through  the  duct. 


(1-3) 


Due  Co  Che  unique  geometry,  Che  corresponding  planes  in  Figure 
1-3  do  noC  separaCe  cells  whose  flow  paccerns  are  idencical,  buc  racher 
are  mirror  images  of  each  ocher.  CompuCacional  difficulcies  or  an  in- 
validacion  of  che  assumpcions  of  periodic  flow  do  noc  resulc,  however. 

The  applicaCion  of  Che  periodic  boundary  condicions  Co  chis  problem  will 
be  discussed  in  Che  Analysis  secdon. 

1.2  Problem  Scacemenc  for  RUFF 

The  usage  of  arcificially  roughened  surfaces  for  enhancing  heac 
transfer  has  been  under  lnvescigacion  for  some  clme  [Bergeles,  1978]. 

This  reporc  will  be  concerned  with  a  type  of  roughness  referred  to  as 
cransverse  or  repeated  ribs. 

The  ribs  are  reccangular  in  shape  and  affixed  to,  or  integral 
with,  a  flat  place.  The  geometry  co  be  investigated  is  similar  to  that 
shown  in  Figure  1-4.  The  dimensions  shown  are  the  rib  height,  h,  the 
rib  length,  1,  the  rib  separation,  s,  the  plate  thickness,  y^,  and  the 
wall  to  centerline  distance,  y^. 

Though  considerable  data  exist  for  repeated  rib  roughness,  there 
remains  uncertainty  in  predicting  the  performance  of  a  particular  re¬ 
peated  rib  surface.  There  have  been  attempts  to  reduce  data  to  a  general 
correlation  [Webb,  1978]  and  to  provide  a  means  for  optimizing  a  rough¬ 
ness  configuration,  [Lewis,  1975  a,b].  These  methods,  however,  provide 
the  designer  with  average  characteristics  of  the  rough  surfaces  and  have 
their  own  uncertainties. 


b)  Dimensioned  Schematic 


Fig.  1-4.  General  Drawings  for  a 
Repeating  Rib  Geometry 


1 


11 


A  heat  exchanger  designer  may  desire  knowledge  of  the  local 
values  of  the  Nusselt  number  and  where  potential  hot  spots  may  occur. 
Additionally,  the  designer  may  wish  to  have  a  convenient  means  of  pre¬ 
diction  and  optimization  at  his  disposal  in  order  to  allow  him  to  survey 
many  designs  without  laborious  usage  of  graphs,  correlations  or  experi¬ 
ments. 

A  general  numerical  program  is  potentially  a  convenient  means  of 
prediction  and  optimization.  The  program  may  be  constructed  to  provide 
the  detailed  information  the  designer  may  desire  as  well.  The  overall 
performance  of  a  repeated-rib  surface  may  be  estimated  by  a  numerical 
solution  of  the  flow  and  heat  transfer  in  a  cell  if  the  entry  region  is 
of  negligible  overall  importance.  Detailed  information  for  much  of  the 
duct  may  be  provided  as  well. 

This  report  continues  an  effort  to  develop  a  numerical  program 
solving  for  the  flow  and  heat  transfer  in  a  cell  begun  by  Short  [1977]. 
The  latest  program  to  be  developed  is  titled  program  RITFF. 

Program  RUFF  utilizes  a  two  cell  approach  as  described  in 
Chapter  3  and  solves  only  for  constant  property,  laminar,  incompressible 
flow.  The  solution  domain  is  shown  in  Figure  1-5.  This  figure  is  ap¬ 
proximately  to  the  scale  of  a  planned  experiment.  Since  the  duct  shown 
in  Figure  1-4  is  symmetric  about  the  centerline,  it  is  only  necessary 
to  solve  from  the  wall  to  the  centerline.  The  energy  equation  is  solved 
both  in  the  fluid  and  the  plate. 
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Fig.  1-5.  Solution  Domain  for  Program  RUFF 
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Extensions  beyond  Short's  work  include  calculation  of  local 
friction  factors  and  Nusselt  numbers,  a  more  flexible  scheme  for  assign¬ 


ing  the  numerical  grid  -  allowing  easier  geometry  definition  -  and  a 
stable  solution  method. 


CHAPTER  2 


ANALYSIS 


2.1  Governing  Equations 

The  problems  considered  in  this  report  are  assumed  to  involve 
steady  state,  recirculating,  laminar  incompressible  flows  with  constant 
properties.  The  fluid  medium  is  presumed  to  be  a  Newtonian  fluid. 

The  governing  equations  for  the  CRACFLO  and  RUFF  program  in  a 
cartesian  coordinate  system  are: 

Continuity: 

3u  .  3v 


(2-1) 


X-momentum: 


f  3u  3v  1  3p  f32u  32u  I  ,  ’it* 

[u3^+v37j’-35  +  w  [§F  +  IF  J  x 


Y-momentum: 
P 

Energy 
pc 


f  3v  . 

3v  1 

3y 

("32v  ,  32v  1 

LU  37  +  v 

3y  J 

w  |jf  +  w  J 

r  3t 

3T  1 

-  k 

P  LU  3x  + 

V  3yJ 

JF  +  sF  J  +  st 

(2-2) 


(2-3) 


(2-4) 


It  has  been  assumed  in  the  energy  equation  that  viscous  dissipation  is 
negligible. 
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The  governing  equations  are  solved  such  that  conditions  of  period¬ 
icity  as  described  by  Short  [1977]  and  separately  by  Patanker,  Liu  and 
Sparrow  [1977]  are  imposed. 

2.2  Equation  Non-Dimensionalization 


Problem  CRACFLO 

The  non-dimensionalization  of  the  momentum  equations  for  problem 
CRACFLO  will  be  illustrated  using  the  x-momentum  equation.  The  non- 


dimensional  variables  are: 


—  x 

X  =  — 

S 


~  _  P_ 

p  puf 


«-3 


B 


—  v 
v  *  — 
U„ 


Introducing  equations  2-5  into  equation  2-2  yields 


(2-5) 


3x"  3 iy  3x  Re 


afu  +  ifu  ]  +  iili 

3 x2  3p  J  PU* 


(2.6) 


where  the  Reynolds  number  is 


Re 


u 


(2.7) 


The  controlling  parameters  in  problem  CRACFLO  may  now  be  seen  to  be  the 
Reynolds  number  and  the  geometry  ratio,  a/g. 
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Problem  RUFF 

The  non-dimensionalization  of  the  momentum  equations  for  problem 
RUFF  will  be  illustrated  using  the  x-momentum  equation.  The  non- 
dimensional  variables  are: 


v  = 


v 


Introducing  equations  2-8  into  equation  2-2  yields 
3x  3y  3x  Re  Zy2  3 y2  ■»  pU| 


(2-9) 


where 


Re 


p0cyll 


(2-10) 


The  energy  equation  is  non-dimensionalized  using  equations  2-8 
and  the  non-dimensional  expression  for  temperature: 

T  -  T 


(2-11) 


qp'yrf 


i 


where  q"  is  Che  imposed  heat  flux,  yfi  the  distance  from  the  plate  wall 

p  ™ 

to  the  centerline  and  kf  the  fluid  thermal  conductivity.  The  non- 
dimensional  energy  equation  in  the  fluid  is 


-  3T  .  -  3T  1 
u  —  +  v  — —  *  - 

3x  3y  Re  Pr 


32T  +  j>^T  +  ST  kf 
aP  sT2  pcpuBq’’ 


(2-12) 


The  non-dimensional  energy  equation  in  the  plate  is  different: 
32T 


Rp/kf  r  -\2~ 


Re  Pr 


a2T 


3y 


c*  *  'V 

bT  f 


(2-13) 


ay2  i  pCptJgq^1 


The  controlling  parameters  that  describe  problem  RUFF  are  seen  to 
be:  the  Reynolds  number,  the  Prandtl  number,  the  duct  geometry  and  the 
ratio,  kp/kf  • 

All  planar  sources  of  heat  such  as  q”  or  q^'  are  non- 
dimensionalized  as  follows : 


q"/q;’ 

Re  Pr 


2.3  Boundary  Conditions  for  CRACFLO 
Two  solution  cells  for  the  CRACFLO  program  are  shown  in  Figure  1-3. 
The  corresponding  planes  are  shown  in  dotted  lines  and  are  numbered  1,  2 
and  3.  The  heavy  solid  lines  represent  impermeable,  stationary  boundaries. 
For  the  solid  boundaries,  the  no-slip  boundary  conditions 

uboundary  "  ® 

(2-14) 

v,  ,  »  0 

boundary 


are  applied. 
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Ic  is  necessary  to  hypothesize  the  form  of  the  velocity  solution 
to  provide  boundary  conditions  at  the  corresponding  planes  1,  2  and  3. 
For  a  two-dimensional  duct,  we  hypothesize  the  velocity  solution  to  be 


u(x,y)  -  u(x,y)encry  + 
region 


u(x,y) 


periodic 


(2-15) 


Since  we  are  interested  only  in  the  solution  for  the  established  vel- 

ocity  field  where  the  effects  of  u(x.y)eni:ry  are  negligible,  ic  is 

region 

u(x,y)  er^0<j^c  which  provides  the  boundary  conditons. 

Periodic  flow  in  this  case  means  steady  flow  with  time,  but  peri¬ 
odic  in  space.  That  is,  the  flow  in  one  solution  cell  is  identical  to 
the  flow  in  all  other  cells  beyond  the  entry  region.  The  corresponding 

planes  separate  the  duct  into  solution  cells.  At  these  corresponding 
planes,  the  hypothesis  of  periodic,  established  flow  results  in  the  follow¬ 
ing  conditions  which  must  be  met  for  the  CRACFLO  problem: 


-v(x,g  +  a)  -  u(g  +  a,  g  +  a  -  x)  =  -v(x  +  g  +  a,0) 

-u(x,g  +  a)  *  v(g  +  a,  g  +  a  -  x)  ■  -u(x  +  g  +  a,0) 

(2-16) 

_3v  (x,g  +  a)  ^  _3v  (g  a,  g  +  a-x)  “  £v  (x  +  g  +■  a>0) 

3y  3x  3y 

3u  (x,g  +  a)  -  3v  (g  +  a,  g  +  a-x)  ■  3u  (x  +-  g  4-  a,0) 

3x  3y  3x 

These  conditions  cannot  be  imposed  directly  because  the  velocities  and 
gradients  are  unknown  and  must  be  determined  as  part  of  the  solution  pro¬ 
cedure. 

The  energy  equation  is  not  solved  in  CRACFLO ,  thereby  removing  the 
need  for  boundary  conditions  on  it. 
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2.4  Boundary  Conditions  for  RUT? 

Two  solution  cells  are  shown, for  program  RUFF  in  Figure  1-5.  The 

heavy  solid  lines  denote  stationary,  impermeable  boundaries  and  the  dot- 

ted  lines  represent  corresponding  planes. 

For  the  solid  boundaries,  the  no  slip  boundary  conditions, 

u  ■  0 

boundary 


(2-17) 


v 


boundary 


0 


are  applied  at  the  fluidplate  interface.  At  the  external  edge  of  the  plate 
(y  »  0)  a  prescribed  heat  flux,  q^'»  imposed. 

At  the  centerline,  we  impose  conditions  of  symmetry: 


3U  (x,y 

3y 


Q. 


0 


p.  +  V  -  0 

3y  r 


£  +  V  •  o 


(2-18) 


To  provide  boundary  conditions  at  the  corresponding  planes,  it  is 
necessary  to  hypothesize  the  form  of  the  velocity  and  temperature  solu¬ 
tion  (Short,  1977].  We  may  hypothesize  solutions  such  as: 


u(*,y)  .  u(*,y)ent.7  -  p,tlodic 

region 

TCx.y)  ■  !(*./> eotry  *  «*  +  T<*.y)p,riodic 

region 


(2-19) 


The  term  mx  in  the  temperature  distribution  accounts  for  the  continuous 
increase  in  temperature  level  due  to  the  imposed  heat  flux. 
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Since  we  are  interested  in  established  flow  where  the  effect  of 
the  entry  region  term  is  negligible,  equations  (2-19)  become: 


u(x,y) ,  , ,  *  u(x,y)  .  , . 

fully  periodic 


T(x,y) 


developed 

fully 

developed 


(2-20) 


mx  +  T(x,y) 


periodic 


The  spatially  periodic  condition  of  the  velocities  implies  that 
flow  in  one  solution  cell  is  identical  to  the  flow  in  all  other  solution 
cells  beyond  the  entrance  region.  The  postulated  periodic  flow  yields 
Che  remaining  boundary  conditions  for  the  two  momentum  equations: 


u(0,y)  » 

u(s  +  w,y)  *  u(2s  + 

2w,y) 

(2-21a) 

v(0,y)  - 

v(s  +  w>y)  *  v(2s  + 

2w,y) 

(2-21b) 

3u  (0,y) 

3u  (s  +w,y)  3u 

(2s  +  2w,y) 

(2-21c) 

3x 

3x 

3x 

3v  (0 ,y) 

3v  (s  +w,y)  m  drv 

(2s  +  2w , y) 

(2-21d) 

3x 

3x 

3x 

The  boundary  conditions  for  the  energy  equation  follow  from  much 
the  same  reasoning,  except  the  addition  of  energy  must  be  treated.  Since 
energy  is  being  added  to  the  fluid,  the  bulk  temperature  of  the  fluid  is 
continuously  rising  so  it  is  impossible  to  relate  scalar  temperatures  in 
the  same  manner  as  scalar  velocities,  equations  (2-21a)  and  (2-21b). 

As  a  result  of  a  constant  heat  flux  applied  to  the  exterior  of  che 
plate,  the  bulk  temperature  of  the  fluid  increases  almost  linearly,  lead¬ 
ing  to  che  term,  mx,  in  the  hypothesized  temperature  solution.  The  con¬ 
ditions  of  periodic  flow  yield  as  a  consequence  identical  temperature 
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profiles  ac  corresponding  planes  when  the  temperatures  are  non- 
dimensionalized  by  che  difference  between  the  bulk  fluid  temperature  and 
the  wall  temperature  at  these  planes 


X  -  T 


T  -  T 

T  -  T 

w 

I  -  T 

„  X_  -  X  | 

B  w 

(2)3  w  \ 

(2-22) 


We  know  that  the  bulk  temperatures  of  the  fluid  are  related  linearly  be¬ 
tween  any  two  corresponding  planes,  i.e.. 


m^s  +  v) 

B©  B® 


T,  *  T_  _  +  +  2w) 

B®  BQ 


(2-23) 


If  we  were  to  assume  chat  I,  and  Xw  are  also  linearly  related  by  the  same 
constant,  m,  equation  (2-22)  is  satisfied  identically.  It  is  che  condi¬ 
tion  of  all  temperatures-  ac  corresponding  planes  being  linearly  related 
that  will  be  used  to  provide  one  streamwise  boundary  condition  on  the 
energy  equation: 

T(0,y)  =  T(s  +  w,y)  -  m(s  +  w)  *  T(2s  +  2w,y)  -  m(2s  +  2w)  (2 -24) 


Another  boundary  condition  at  the  corresponding  planes  is  pro¬ 
vided  by  differentiating  equation  (2-24)  with  respect  to  x 


^.(0,y)  a  +  (2s  +  2w,y) 

3x  3x  3x 


(2-25) 


CHAPTER  3 


NUMERICAL  TECHNIQUE 

The  primary  difficulty  in  solving  problems  with  periodic  boun¬ 
dary  conditions  is  that  the  values  of  velocity  and  temperature  on  the 
upstream  and  downstream  corresponding  planes  are  unknown.  Typically 
all  values  and/or  fluxes  are  known  on  the  boundaries  in  numerical  prob¬ 
lems.  This  chapter  concerns  itself  with  the  method  used  in  programs 
CRACFLO  and  RUFF  to  impose  periodic  boundary  conditions. 

3.1  History 

There  are  a  number  of  numerical  techniques  to  impose  periodic 
boundary  conditions.  All  methods  have  a  common  action.  The  methods 
act  upon  those  boundary  velocities  which  are  subject  to  periodic  boun¬ 
dary  conditions  in  an  orderly  fashion  seeking  those  boundary  values 
which  are  correct  according  to  the  physics  of  the  problem.  The  boundary 
temperatures  are  similarly  treated  in  cases  where  the  energy  equation  is 
solved. 

A  total  of  five  methods  have  been  investigated  by  the  author. 
These  will  be  briefly  described  in  the  chronological  order  of  their  in¬ 
vestigation.  The  term  "updating,"  used  in  the  following  discussion  is 
defined  as  any  rational  process  by  which  the  upstream  and  downstream 
boundary  conditions  are  altered  in  an  attempt  to  affect  a  solution. 

The  first  technique  was  the  one-cell  successive  substitution 
method  of  Short  [1977],  This  method  was  found  to  be  unstable  in  the 
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sense  chat  it  produced  a  divergent  solution.  Later,  this  instability 
was  determined  to  be  due  to  updating  the  upstream  and  downstream  boun¬ 
dary  conditions  every  iteration. 

The  second  technique  was  a  sensitivity  matrix  method  invented 
by  the  author  but  related  to  Newton’s  method  [Cebeci  and  Bradshaw, 

1977,  sec.  4.2.2].  This  method  utilized  a  matrix  of  sensitivity  co¬ 
efficients  produced  by  perturbing  a  test  flow  field  for  a  particular 
geometry.  The  sensitivity  matrix  provided  linear  relationships  be¬ 
tween  changes  in  the  boundary  velocity  values  and  changes  in  the  vel¬ 
ocities  on  the  internal  plane.  These  relationships  furnished  a  means 
of  updating  to  obtain  periodic  boundary  conditions.  The  sensitivity 
matrix  method  required  more  accuracy  and  execution  time  of  the  com¬ 
puter  program  than  appeared  economically  feasible. 

The  third  technique  was  the  internal  plane  method.  This  method 
utilizes  the  velocities  on  the  internal  plane  directly  to  provide  new 
values  on  the  boundary  planes.  This  method  is  stable  and  applicable 
to  both  the  CRACFLO  and  RUFF  problems.  One  drawback  is  the  require¬ 
ment  for  considerable  execution  time.  However,  this  method  was  sel¬ 
ected  to  solve  the  CRACFLO  and  RUFF  problems,  partly  due  to  its  reli¬ 
ability. 

The  fourth  technique  is  called  a  stabilized  one-cell  succes¬ 
sive  substitution  method.  This  is  the  method  of  Short  [1977],  stabil¬ 
ized  by  updating  only  every  ten  iterations.  Though  the  solution 
domain  is  half  that  of  the  internal  plane  update  method,  execution 
time  Is  greater  due  to  a  lower  convergence  rate. 


The  fifth  technique  is  that  of  Patankar,  Liu  and  Sparrow  [1977], 
This  method  utilizes  a  cyclic  tridiagonal  matrix  in  the  streamwise  di¬ 
rection  to  impose  periodic  boundary  conditions  directly  as  part  of  the 
solution  of  the  finite  difference  equations  rather  than  separately  as  is 
done  in  the  first  four  methods.  The  author  was  unable  to  construct  an 
operable  version  of  this  method  for  testing.  While  the  application  of 
this  method  to  CRACFLO  presents  difficulties,  the  cyclic  tridiagonal 
matrix  method  may  be  an  excellent  alternative  for  RUFF.  However,  since 
the  TEACH  program  uses  an  upwind  differencing  procedure  for  strong  con¬ 
vection,  this  technique  may  not  be  effective  when  there  is  dominant  out¬ 
flow  at  the  downstream  boundary.  The  Internal  Plane  Method,  abbreviated 
IPM,  will  now  be  discussed  in  detail.  The  other  four  methods  are  cover¬ 
ed  in  more  detail  by  Faas  [1979;  Appendix  D]. 

3.2  Internal  Plane  Method 

The  IP  method  Is,  in  a  sense,  a  fusion  of  the  method  of  Short 
[1977]  and  certain  elements  of  the  sensitivity  matrix  method  mentioned 
above.  The  IP  method  uses  two  cells  sharing  a  common  corresponding 
plane  and  requires  the  external  corresponding  planes  to  have  identical 
scalar  values  of  velocity.  It  also  uses  the  internal  plane  to  provide 
information  leading  to  the  estimation  of  new  boundary  values. 

The  IP  method  will  be  treated  in  detail  using  the  solution  do¬ 
main  of  program  CRACFLO  as  an  example.  The  geometry  is  repeated  in 
Fig.  3-1.  The  dotted  lines  numbered  1,  2  and  3  are  the  corresponding 
planes.  Planes  1  and  3  are  the  boundary  planes  and  plane  2  is  the  inter¬ 
nal  plane.  The  internal  plane  is  a  mirror  image  of  the  boundary  planes, 
with  point  A  corresponding  to  point  A'. 
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Fig.  3-1.  CRACFLO  Geometry 

The  general  numerical  technique  of  IPM  is  to  estimate  initial  val¬ 
ues  of  u  and  v  velocities  at  the  planes  1  and  3  and  to  solve  for  the  flow 
field  numerically.  The  influence  of  the  duct  geometry  as  well  as  the 
boundary  conditions  determines  the  u  and  v  velocities  at  plane  2.  Hence, 
these  u  and  v  velocity  values  become  better  estimates  of  the  actual  vel¬ 
ocity  values  that  satisfy  the  periodic  boundary  conditions. 

These  deduced  values  are  then  substituted  for  the  previous  boun¬ 
dary  velocity  values  and  the  process  is  repeated  until  periodic  boundary 
conditions  are  met.  The  process  of  starting  with  a  set  of  boundary  values 
and  ending  with  their  replacement  with  values  from  the  internal  plane  will 
be  referred  to  as  a  "cycle."  This  process  is  distinct  from  an  iteration 
of  the  core  program. 

Many  details  have  been  omitted  above  and  require  further  discus¬ 
sion.  A  simple  model  using  four  nodes  cross-stream  will  be  used  for  pur¬ 
poses  of  illustration.  Figure  3-2  shows  the  location  and  title  of  the 
various  pertinent  velocities.  The  subscript  "b"  denotes  a  boundary  vel¬ 
ocity,  the  subscript  "int"  an  internal  plane  velocity,  the  superscript 


ocity  locations  and  titles  for 
llustrative  example 


"a"  refers  to  an  upstream  boundary,  and  the  superscript  "b"  to  a  down¬ 
stream  velocity.  The  superscripts  "a"  and  "b"  are  necessary  due  to  the 
staggered  grid  utilized  by  TEACH,  the  program  which  is  a  basis  of 
CRACFLO  and  RUFF. 

The  first  step  taken  by  program  CRACFLO  is  to  initialize  the 
variable  fields  and  the  boundary  conditions.  The  u  and  v  velocities  are 
initially  set  to  zero  except  for  the  boundary  velocities  ^  which  are 
initialized  to  represent  a  parabolic  velocity  profile  whose  flow  direc¬ 
tion  is  the  same  as  shown  in  Figure  3-2.  The  subscript  "n"  may  take  on 
the  values  1,  2,  3  or  4  in  this  example. 

The  core  program  uses  the  initial  boundary  conditions  to  solve 
for  u  and  v  velocities  and  pressures  within  the  solution  domain  until  a 
convergence  criterion  ym  is  met. 

The  core  program  utilizes  a  convergence  criterion  based  on  the 
calculated  mass  flow  imbalance  over  the  solution  domain.  A  mass  flow 
imbalance  of  a  nodal  control  volume  is  referred  to  as  its  residual  mass 
source.  The  core  program  is  considered  converged  when  the  absolute  sum 
of  the  mass  residual  sources  of  the  nodal  control  volumes  is  less  than 
a  certain  fraction  of  the  mass  flow  rate  through  the  solution  domain. 
Mathematically  this  is 
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where  ym  is  the  chosen  fraction  of  the  mass  flow  rate. 
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Core  program  convergence  is  complicated  by  the  ability  to  achieve 
a  desired  tolerance  value  of  Ym  when  the  estimated  boundary  velocity  val¬ 
ues  are  far  from  those  of  periodic  flow.  This  situation  is  almost  always 
the  case  for  the  initial  boundary  values.  This  complication  is  solved  by 
selecting  an  initial  tolerance,  Ym^  ^  ,  which  is  sufficiently  large  to 
allow  convergence  to  be  achieved  with  initial  boundary  values.  This  value 
is  reduced  as  the  program  proceeds  by  a  preset  constant,  R,  prior  to  the 
commencement  of  the  next  cycle.  The  tolerance  is  reduced  until  a  mini¬ 
mum  tolerance,  Ym  .  ,  is  reached  or  exceeded.  The  value  of  Ym  is  then 
min 

fixed  at  Ym^^  for  all  subsequent  cycles . 

Once  the  core  program  has  converged,  the  program  then  checks  for 
overall  convergence.  Overall  convergence  is  determined  by  comparing  how 
well  the  deduced  boundary  values  satisfy  the  periodic  boundary  conditions. 
Ideally,  one  condition  to  be  met  is 
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int 


(3-2) 


where  the  subscript  "b"  refers  to  a  value  on  the  boundary  planes  and  the 
subscript  "int"  refers  to  a  value  on  the  internal  plane.  Since  the  in¬ 
ternal  plane  is  rotated  90°  with  respect  to  the  boundary  plane,  we  must 
compare  u  velocities  on  one  plane  with  v  velocities  on  the  other  with 
the  appropriate  correction  of  sign. 

Experience  with  the  program  shows  that  the  streamwise  velocities 
apparently  dominate  the  physics  of  the  flow  to  such  an  extent  that  sat¬ 
isfying  equation  (3-2)  results  in  the  satisfaction  of  the  remaining  peri¬ 
odic  boundary  conditions  as  detailed  in  chapter  two.  This  dominance  is 
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exploited  co  provide  overall  convergence  criteria  which  are  relatively 
simple  to  apply.  Overall  convergence  is  considered  to  be  achieved  when 
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where  the  subscript  "a"  refers  to  a  particular  corresponding  node  on  the 
internal  and  boundary  planes  and  Yv  is  a  tolerence  value.  This  criterion 
must  be  met  by  every  pair  of  corresponding  velocities  before  overall  con¬ 
vergence  is  considered  accomplished.  If  overall  convergence  is  accom¬ 
plished,  then  the  program  prints  the  variable  arrays  and  local  flow  angles 
then  stops.  If  the  convergence  criteria  are  not  met,  the  program  proceeds 
to  perform  an  update  of  the  trial  boundary  values. 

The  program  numerically  integrates  the  velocity  profile  at  the 

internal  plane  to  determine  the  mass  flow  rate,  m.  ,  across  it.  A  flow 
^  mt 

correction  factor,  F  ,  is  calculated  bv 
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where  ®r.e£  is  calculated  from  the  Reynolds  number  and  fluid  properties. 

Next  the  velocities  at  the  internal  plane  are  used  to  replace 

chose  at  the  boundaries  as  follows: 
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where  it  may  be  seen  the  purpose  of  Fc  is  to  scale  the  velocity  profile 
at  the  internal  plane  to  conserve  the  overall  mass  flow  rate. 

The  correction  factor,  F  ,  is  printed  every  update  as  an  indica¬ 
tion  of  how  the  solution  procedure  is  working.  The  correction  factor 
should  approach  a  constant  value  of  one  as  the  program  moves  toward 
overall  convergence.  This  is  observed  in  all  cases. 

The  program  then  reduces  the  value  of  the  core  program  convergence 

criteria,  y  ,  by  a  factor  of  R.  If  this  results  in  the  value  of  y  being 
m  tn 

reduced  below  a  minimum  value,  Ym  .  ,  the  value  of  Y  is  set  to  Tm  ,  . 

min  m  min 

The  core  program  then  solves  for  the  u  and  v  velocities  and  pres¬ 
sures  with  the  new  boundary  values  on  the  boundary  planes.  The  process 
is  repeated  until  overall  convergence  is  reached.  A  flow  chart  of  the 
entire  process  is  shown  in  Figure  3-3. 


3.3  Application  of  the  IP  Method  to  RUFF 
The  IP  technique  is  applied  to  RUFF  in  a  straight-forward  fashion 
comparable  to  CRACFLO  insofar  as  the  velocities  are  concerned.  Additional 
calculations  are  required  to  enforce  periodic  boundary  conditions  on  the 
energy  equation  as  detailed  in  chapter  two. 


Overall  convergence  of  RUFF  is  reached  when  the  following  is  at¬ 


tained 
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In  order  to  maintain  a  constant  mass  flow  rate,  the  streanwise 
velocities  are  corrected  at  each  update  by  a  correction  factor  Fc  as 
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defined  in  equation  (3-4).  The  update  performed  is  shown  below 


u  ■  F  u. 
b,n  c  mt,n 

a  _  a 
vb,n  vint,n 

b  b 

b,n  int,n 


(3-9) 

(3-10) 

(3-11) 


The  superscripts  "a”  and  "b"  refer  to  the  upstream  and  downstream 
boundary  and  internal  values  respectively  as  related  in  the  CRACFLO  ex¬ 
ample  and  shown  in  Figure  3-4. 

A  constant  heat  flux  is  applied  to  the  outer  edge  of  the  ribbed 
plate  in  program  RUFF.  This  requires  the  bulk  temperature  of  the  fluid 
between  corresponding  planes  to  increase  with  x  by  the  incremental  amount, 
m,  defined  by 
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PVW 


(3-12) 


where  q^  is  the  imposed  heat  flux,  c^  the  fluid  specific  heat,  Tg  the 
bulk  temperature  of  the  fluid  and  U_  the  bulk  fluid  velocity. 

D 

The  specification  of  the  fluid  bulk  temperature  at  any  section 
of  the  computational  cell  means  the  bulk  temperature  at  any  corresponding 
plane  of  the  duct  is  also  specified.  Selecting  an  inlet  bulk  temperature 
for  an  upstream  boundary  condition  automatically  determines  the  bulk  tem¬ 
perature  at  the  downstream  boundary.  These  bulk  temperatures  must  remain 
fixed  regardless  of  the  boundary  temperature  profiles  to  maintain  overall 
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Fig.  3-4.  Location  of  the  "a"  and  "b"  Planes 
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conservation  of  energy.  The  boundary  temperature  profiles  are  updated 
such  that  the  inlet  and  outlet  bulk  temperatures  conserve  energy  as  fol¬ 
lows. 

When  Ruff  initializes  the  temperature  field  for  a  particular  prob¬ 
lem,  the  bulk  temperatures  at  the  two  internal  corresponding  planes,  T^a 
and  T_,  ,  are  recorded.  Two  internal  planes  are  required  for  the  tempera- 

DO 

tures  as  was  necessary  for  the  v-velocities  because  of  the  staggered  grid 
used  in  program  TEACH. 

The  energy  equation  is  solved  in  each  cycle  after  the  momentum 
equations  are  solved  for  the  velocity  field.  This  is  possible  since  the 
fluid  properties  are  not  temperature  dependent  in  program  RUFF;  thereby 
the  momentum  equations  are  decoupled  from  the  energy  equation.  The  energy 
equation  is  solved  to  an  accuracy  where  the  sum  of  the  individual  nodal 
energy  imbalances  is  less  than  a  specified  fraction  of  the  imposed  heat 
transfer  through  the  wall,  that  is. 
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where  is  the  tolerance  fraction  and  &x  -  2(S  +  W) .  (In  the  present 
version  of  the  program  the  units  are  incorrect  on  the  righthand  side  due 
to  an  oversight  when  constructing  the  code.  Operationally,  this  diffi¬ 
culty  has  had  little  effect  on  program  performance.) 

The  temperature  profiles  at  the  internal  corresponding  planes 
are  integrated  to  determine  the  bulk  temperatures,  T0^  and  T^*  at  those 
planes.  Temperature  correction  factors,  T^  and  T^,  are  then  produced 
using  the  bulk  temperatures,  T0a  and  T^,  deduced  when  RUFF  was  initial¬ 


ized: 
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The  temperature  boundary  conditions  are  then 

updated  upstream  and 

downstream  by  the  following  process: 

a  rp 
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A  general  flowchart  for  RUFF  is  shown  in  Figure  3-5. 


3.4  Confidence  Tests 

Program  CRACFLO 

Several  tests  of  program  CRACFLO  were  conducted  to  determine  a) 
the  uniqueness  of  the  numerical  solutions  found,  and  b)  the  effect  of  re¬ 
fining  the  grid  mesh.  The  latter  test  is  also  known  as  a  test  of  grid 
independence. 

When  discussing  the  tests,  a  convenient  way  of  expressing  the 
grid  refinement  is  by  the  number  of  nodes  cross-stream,  The  number 

of  nodes  cross-stream  is  the  number  of  nodes  on  a  corresponding  plane  that 
lie  between  the  solid  walls.  For  example,  the  number  of  nodes  cross¬ 
stream  shown  in  Figure  3-2  is  four.  (When  a  particular  number  of  nodes 
cross-stream  is  mentioned  in  the  following,  the  reference  is  to  a  specif¬ 
ic  grid  scheme  that  was  used  and  not  to  one  of  many  other  grid  schemes 


that  satisfy  the  requirement  of  having  the  specified  number  of  nodes 
cross-stream.) 

The  uniqueness  tests  were  conducted  with  a  grid  having  12  nodes 
cross-stream,  a  Reynolds  number  of  100  and  a  geometry  ratio,  a/y  *  0.5. 
Instead  of  the  usual  parabolic  profile  as  initial  inlet  and  outlet  bound¬ 
ary  velocity  profiles,  a  sine  curve  and  slug  flow  profile  were  used.  The 
direction  of  the  flow  was  identical  in  all  three  cases.  All  three  cases 
yielded  identical  flow  fields,  pressures  and  average  friction  factors. 

Additionally,  a  slug  profile  with  reversed  sign  was  used  at  the 
inlet  and  outlet  as  an  initial  velocity  profile.  That  is,  the  direction 
of  the  flow  through  the  solution  domain  was  reversed.  The  Reynolds  num¬ 
ber  was  identical  with  the  previous  tests.  The  resulting  flow  field  was 
the  inverse  of  the  previous  flow  fields  mentioned  above  and  the  average 
friction  factor  was  identical.  Program  CRACFLO,  to  the  extent  it  was 
tested,  exhibited  the  ability  to  attain  a  unique  solution  for  a  unique 
Reynolds  number  and  geometry  ratio. 

Solutions  were  found  for  Re  *  100,  a/g  *  0.5  and  6,  12,  24  and  36 
nodes  cross-stream  for  a  grid  independence  test.  The  solution  for  36 
nodes  cross-stream  did  not  reach  full  convergence;  it  had  progressed  far 
enough  to  determine  the  average  non-dimensional  pressure  drop  between 
corresponding  planes  but  not  the  finest  details  of  the  flow. 

The  values  of  the  average  non-dimensional  pressure  drop  between 
corresponding  planes  versus  the  number  of  nodes  cross-stream  for  the  above 
solutions  are  plotted  in  Figure  3-6.  They  show  a  definite  trend  towards 
an  asymptotic  value  with  increasing  grid  refinement,  indicating  a  grid 


dimensional  Pressure  Drop  Between  Corresponding  Planes 
as  a  Function  of  the  Number  of  Nodes  Cross-stream 
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independent  solution  can  be  obtained.  This  manner  of  plotting  results 
of  the  grid  independence  demands  more  of  the  numerical  technique  in  per¬ 
formance  than  other  means  of  plotting  the  solution  results  such  as  against 
the  total  number  of  nodes  in  the  solution  domain.  Hence  the  trend  indi¬ 
cated  in  the  plot  is  stronger  evidence  than  may  first  appear. 

Plotted  as  well  in  Figure  3-5  are  the  results  for  a  grid  indepen¬ 
dence  test  for  Re  =  10  and  a/ g  =  0.5.  Results  of  both  grid  independence 
tests  are  summarized  in  Table  3-1  and  Table  3-2.  Additionally,  velocity 
profiles  from  the  internal  corresponding  planes  for  the  solution  for 
Re  31  100,  a/g  *  0.5,  are  plotted  in  Figure  3-7.  Due  to  the  staggered 
grid  used  in  program  TEACH,  there  are  no  values  plotted  in  Figure  3-6 
for  y/g  =0.0  and  y/g  =  1.0.  Any  points  apparently  at  these  surfaces 
are  actually  for  a  y/g  value  very  close  to  the  wall.  To  emphasize  the 
region  near  the  wall,  a  non-uniform  grid  was  used  in  all  cases. 

The  results  of  the  grid  independence  tests  indicate  that  a  grid 
having  12  nodes  cross-stream  is  sufficient  to  produce  accurate  average 
non-dimensional  pressure  drops  for  Re  10.  For  the  results  presented, 
a  grid  having  24  nodes  cross-stream  is  used  for  Re  >  10  as  a  compromise 
between  accuracy  and  economy. 

The  necessity  for  using  more  nodes  cross-stream  for  Re  >  10  seems 
to  stem  from  the  growth  of  the  recirculation  zones.  Above  Re  =  10,  the 
recirculation/separation  zones  grow  to  intersect  and  cross  the  correspond¬ 
ing  planes.  The  appearance  of  this  reversed  flow  at  these  planes  injects 
a  flow  structure  of  finer  scale  than  a  12  nodes  cross-stream  grid  can 
adequately  model.  The  solution  is  to  refine  the  grid,  in  this  case  to  a 
grid  of  24  nodes  cross-stream. 


Velocity  Profiles  on  the  Internal  Corresponding  Plane  for  Grid  Independence 

Tests,  Re  »  100  and  a/g  *»  0.5 
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Program  RUFF 

A  limited  test  of  program  RUFF  was  performed.  This  test  was  to 
execute  program  RUFF  with  the  rib  height  set  to  zero,  i.e.,  as  a  parallel 
plate  geometry.  Analytical  values  of  the  Nusselt  number  and  friction 
factor  for  laminar  flow  between  parallel  plates  are  widely  available  in 
existing  literature  [Kays,  1966].  These  are: 


Nu  =  8.235 
6 


(3-18) 


f  = 


Re 


where  the  Reynolds  number  in  equation  (3-18)  is  based  on  half  the  distance 


between  the  plates  (one  fourth  of  the  hydraulic  diameter)  as  it  is  in  pro¬ 
gram  RUFF. 

Program  RUFF  was  executed  with  Re  =  100  and  a  grid  of  40  x  25 
nodes,  of  which  40  x  18  are  in  the  fluid  and  the  remainder  are  in  the 
solid  plate.  The  results  were: 


Nu  =  8.28 
f  =  0.0604 

Nu  “8.28 
w 

f  =0.06 
w 

where  the  subscript  w  indicates  a  local  value  at  the  wall.  Only  one  local 
value  of  each  is  tabulated;  all  were  essentially  identical. 

The  value  of  the  local  and  average  Nusselt  number  is  0.5%  larger 
than  the  analytical  value.  This  is  considered  acceptable  by  the  author. 

The  friction  factor,  both  local  and  average,  calculated  from 


equation  (3-18)  is  0.06.  The  local  friction  factors  calculated  by  pro¬ 
gram  RUFF  agree  very  well,  varying  only  in  the  fifth  significant  figure. 
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The  average  friction  factor  has  a  similar  accuracy.  The  error  for  both 
local  and  average  friction  factors  is  0.7%  and  is  considered  acceptable 
by  the  author. 

A  solution  for  laminar  flow  in  a  duct  with  repeated  rib  roughness 
was  generated  following  this  test.  The  results  are  presented  in  the  next 
chapter.  The  grid  used  had  non-uniform  spacing  with  40  nodes  in  the 
x-direction  and  25  nodes  in  the  y-direction.  The  region  between  the  plate 
and  centerline  was  allocated  18  nodes,  in  the  y-direction,  with  5  nodes  al¬ 
located  over  the  rib  height. 


3.5  Stability  and  Convergence 

Neither  program  exhibited  any  tendency  towards  instability  under 
any  conditions  used  in  the  course  of  the  author's  research.  The  core 
program  has  its  own  stability  requirement  as  described  by  Faas  [1979; 
Appendix  A]  and  was  stable  under  all  conditions. 

Program  CRACFLO  exhibited  the  ability  to  converge  towards  a  solu¬ 
tion.  The  convergence  rate,  however,  was  lower  than  the  author  preferred. 
The  execution  time  with  a  grid  of  24  nodes  cross-stream,  Re  =  100  and 
a/g  ■  0.5  was  188  seconds  or  0.11  seconds  per  node  on  a  CDC  CYBER  175. 

The  execution  time  with  12  nodes  cross-stream.  Re  =  10,  and  a/g  =  0.5  was 
27.7  seconds  or  0.047  seconds  per  node. 

Program  RUFF  also  exhibited  the  ability  to  converge  uniformly  to 
a  solution,  but  had  an  even  slower  convergence  rate.  The  execution  time 
for  the  problem  investigated  was  257  seconds  or  0.29  seconds  per  node  on 
a  CDC  CYBER  175.  The  24  nodes  cross-stream  case  for  CRACFLO  mentioned 
above  had  a  grid  of  57  x  33  compared  to  RUFF  which  has  40  x  25. 


43 


3.6  The  Core  Program 

The  core  program  used  in  both  program  CRACFLO  and  RUFF  is  program 
TEACH  developed  by  Gosraan  at  the  Imperial  College,  London,  England.  Pro¬ 
gram  TEACH  is  documented  by  Gosman  and  Ideriah  [1976]  and  Gosman  and  Pun 
[1973]  and  is  described  in  the  open  literature  by  Khalil,  Spalding  and 
Whitelaw  [1975],  as  well  as  being  extensively  documented  by  Faas  [1979]. 

In  general,  TEACH  solves  the  momentum  and  energy  equations  in 
terms  of  their  primative  variables  for  two  dimensional,  recirculating 
flows.  It  is  necessary  to  have  a  program  capable  of  handling  recirculat¬ 
ing  flows  as  they  often  occur  for  rough  walls.  A  further  advantage  of 
program  TEACH  is  believed  to  be  its  usage  of  a  hybrid  differencing  tech¬ 
nique  which  uses  central  or  upwind  differencing  locally  according  to  the 
local  flow  conditions  resulting  in  greater  overall  accuracy. 

Modifications  of  program  TEACH  for  the  author's  research  were 
minimal.  The  only  change  of  note  was  the  addition  of  FORTRAN  code  to 
subroutine  LISOLV.  This  code  allows  LISOLV  to  perform  alternating  line- 
by-line  solutions  sweeping  from  left  to  right  and  then  top  to  bottom. 

The  change  was  originally  made  in  program  CRACFLO  because  the  main 
flow  direction  of  the  fluid  was  from  left  to  right  and  top  to  bottom, 
but  was  retained  for  RUFF  when  it  was  found  it  increased  the  convergence 


rate. 


CHAPTER  4 


RESULTS  AND  DISCUSSION 

4. 1  Problem  CRACFLO 

The  Reynolds  number  based  on  Che  gap  width  g  and  the  geometry 
ratio,  a/g,  are  the  controlling  parameters  for  problem  CRACFLO.  Utilizing 
the  IP  method,  we  constructed  program  CRACFLO  and  flow  solutions  were 
produced  wherein  these  parameters  were  varied.  The  average  friction 
factor  based  on  the  pressure  drop  between  corresponding  planes  was  deter¬ 
mined  for  each  pair  of  Reynolds  number  and  geometry  ratio  investigated. 

Definitions 

The  average  friction  factor  used  in  this  report  is  defined  as 

AP 

£--So7Tf  (4-1) 


where  g  is  the  gap  width,  L^  the  mean  distance  between  corresponding 
planes  and  AP  is  the  average  pressure  drop  between  corresponding 
planes  defined  by 


AP 

cp 


(4-2) 


where  the  circled  number  subscript  indicates  the  corresponding  plane 
where  the  integration  is  performed  and  l  is  a  length  variable  aligned 
with  the  plane.  In  the  program  the  pressure  drop  between  the  next  two 
corresponding  surfaces  is  also  calculated  to  provide  insight  into  the 
progress  of  the  solution  procedure. 
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The  distance  L  is  defined  in  Figure  4-1  as  the  distance  along  the 
m 

dashed  line  connecting  the  intersecting  centerplanes  from  A  to  A' .  For 
the  case  illustrated,  in  Figure  4-1: 


Fig.  4-1.  Definition  of  L 

m 


One  may  note  that  as  a  -*00,  the  effect  on  the  average  friction  factor 
caused  by  the  comer  regions,  which  correspond  to  the  constant  1.0  in 
equation  (4-3),  becomes  smaller  and  smaller.  Hence  as  a  approaches00, 
the  average  friction  factor  should  approach  that  for  laminar  flow  between 
parallel  planes.  The  lower  limit  for  a  in  the  present  version  of  CRACFLO 
is  a  value  of  0.0. 

Several  geometries  for  a  corrugated  channel  are  presented  in 
Figure  4-2  to  give  the  reader  an  impression  of  the  various  shapes  which 
can  be  calculated  by  the  program.  As  a/g  increases,  the  situation  ap¬ 
proaches  a  series  of  parallel  plate  ducts  connected  at  right-angle  bends. 

The  Reynolds  number  used  in  this  report  is  defined  by: 


Re  - 


u 


(4-4) 
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This  definition  was  chosen  because  g  was  selected  as  the  characteristic 
length  by  which  all  other  lengths  would  be  non-dimensionalized  in  problem 
C RAC FLO . 

Effect  of  Reynolds  number  on  f  with  a/g  held  constant 

A  series  of  computer  solutions  yielding  average  friction  factors 
was  performed  holding  the  geometry  ratio  constant  while  varying  the  Reyn¬ 
olds  number.  Two  geometries  were  investigated;  one  having  a  geometry 
ratio  of  0.5  and  the  other  a  geometry  ratio  of  0.0.  The  latter  repre¬ 
sents  the  lower  limit  of  geometry  ratios  which  can  be  calculated. 

The  average  friction  factors  for  these  two  geometries  are  plotted 
against  the  Reynolds  number  in  Figures  4-3  and  4-4.  Noticeable  in  both 
figures  is  the  non-linear  shape  of  the  curves  (on  logarithmic  coordinates) 
and  the  presence  of  a  minimum.  The  typical  behavior  of  laminar  flows  in 
ducts  yields  friction  factor  curves  that  are  inversely  proportional  to  the 
Reynolds  number  so  they  form  straight  lines  with  negative  slopes  on  log¬ 
arithmic  coordinates  as  shown  in  Figure  4-4  [Knudsen  and  Katz,  1958]. 

Laminar  wake  flow  about  a  circular  cylinder  -  as  well  as  other 
blunt  bodies  -  has  a  drag  coefficient  curve,  shown  in  Figure  4-6,  similar 
to  the  curve  in  Figure  4-2.  (However,  Figure  4-6  is  a  semi-logarithmic 
plot.)  Figure  4-5  also  relates  the  behavior  of  the  coefficient  of  drag 
to  the  recirculation  zones  or  wake  behind  the  cylinder.  This  observation 
suggests  the  possibility  that  recirculation  zones  in  the  corrugated  duct 
may  control  the  behavior  of  the  average  friction  factor  at  higher  Reynolds 


numbers . 


Fig.  4-3.  Average  Friction  Factor  Variation 
for  a  Geometry  Ratio  of  0.5 
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Re 


Average  Friction  Factor  Variation 
for  a  Geometry  Ratio  of  0.0 


Fig.  4-4. 


num*>e 


Fig.  4-5.  Laminar  Flow  in  Circular  Tubes 
[Knudsen  and  Katz,  1958] 


Fig.  4-6.  Flow  About  a  Circular 
Cylinder  [Murdock,  1976] 


The  effect  of  the  recirculation  zones 

The  effect  of  the  recirculation  zones  will  be  investigated  at  a 
constant  geometry  ratio,  a/g  *  0.5.  A  series  of  diagrams  indicating  the 
local  flow  angles  in  the  corrugated  channel  aids  in  interpreting  the  size 
and  effect  of  the  recirculation  zones.  These  flow  angle  diagrams  are 
shown  in  Figures  4-2,  4-8,  4-9  and  4-10.  All  flow  diagrams  are  plotted 
from  the  results  of  a  12-node  cross-stream  grid.  A  finer  grid  solution 
is  not  warranted  as  the  gross  details  are  of  primary  interest  here. 

The  flow  field  in  Figure  4-6  is  representative  of  the  flow  for 
Re  <_  10.  The  recirculation  zones  are  almost  non-existant ,  hence  the  flow 
follows  the  boundaries  closely  and  moves  uniformly  downstream  with  little 
or  no  separation  or  reversed  flow.  This  flow  structure  is  similar  to 
Stokes'  flow  and  it  is,  therefore,  not  surprising  that  the  average  fric¬ 
tion  factor  curve  is  nearly  linear  for  Re  <  10  on  a  logarithmic  plot 
[Schlichting,  1968]. 

The  recirculation  zones  grow  as  the  Reynolds  number  is  increased. 
At  Re  =  25  the  friction  factor  curve  has  diverged  from  its  linear  appear¬ 
ance.  The  recirculation  zones  have  grown  to  over  1/2  z  in  diameter  as 
shown  in  Figure  4-8.  A  small  secondary  recirculation  zone,  represented  by 
a  single  arrow,  has  appeared  between  the  major  eddy  and  the  corner.  A 
small  separation  region  also  occurs  behind  the  convex  corner  over  which 
the  fluid  flows. 

Figure  4-9  shows  the  flow  field  concurrent  or  nearly  concurrent 
with  the  minimum  in  the  average  friction  factor.  The  primary  recircula¬ 
tion  zones  have  grown  to  a  diameter  of  approximately  3/4  g  and  are  con¬ 
nected  with  the  separated  flow  behind  the  convex  corner.  The  secondary 
recirculation  zones  have  grown  as  well  but  are  still  small. 
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A  flow  field  for  a  Reynolds  number  above  that  corresponding  to 
the  minimum  average  friction  factor  is  shown  in  Figure  4-10.  The  primary 
recirculation  zones  have  grown  slightly  and  little  change  is  observable 
in  the  small  secondary  recirculation  zones. 

The  flow  field  for  Re  >  100  is  represented  by  Figure  4-11.  The 
primary  recirculation  zones  have  not  changed  noticeably,  but  the  secon¬ 
dary  recirculation  zone  has  grown  about  150  percent  over  its  size  in 
Figure  4-9.  In  both  cases  the  main  through  flow  streamlines  appear  to 
cross  the  duct  relatively  directly  and  to  concentrate  near  the  far  wall 
after  passing  over  a  comer. 

A  relationship  between  the  average  friction  factor  and  the  recir¬ 
culation  zones  may  now  be  inferred.  The  recirculation  zones  are  small  or 
non-existent  at  low  Reynolds  numbers,  consequently  the  average  friction 
behaves  like  the  friction  factor  for  Stokes'  flow.  The  recirculation 
zones  grow  as  the  Reynolds  number  increases  and  begin  to  affect  the  over¬ 
all  flow.  This  effect  is  characterized  by  the  decrease  in  slope  of  the 
average  friction  factor  curve. 

The  average  friction  factor  reaches  a  minimum  approximately  con¬ 
current  with  the  stabilization  of  the  primary  recirculation  zone  size. 
Hereafter,  the  secondary  recirculation  zone  grown  preferentially  and 
appears  to  account  for  the  behavior  of  the  average  friction  factor  above 
Re  %  100.  At  the  higher  Reynolds  numbers  the  main  throughflow  becomes 
more  concentrated  and  resembles  a  two-dimensional  jet  passing  through  a 
series  of  expansions  and  contractions  or  connected  asymmetric  2-D  nozzles 
and  diffusers.  The  approximately  constant  friction  factor  corresponds  to 
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the  loss  coefficient  for  comparable  geometries  at  high  Reynolds  numbers 
[Vennard  and  Street,  1975,  sec.  9.14]. 

The  presence  of  these  recirculation  zones  may  account  for  the 
differences  between  Figure  4-3  and  the  upper  curve  of  Figure  4-12  excerpt¬ 
ed  from  Kays  and  London  [1964].  This  w,  vy  fin  geometry  was  selected 
since  it  most  closely  resembles  the  geometry  which  produced  Figure  4-3, 
i.e.,  a  geometry  ratio  of  0.5.  The  fins  of  Figure  4-12  would  corres¬ 
pond  to  a  geometry  ratio  of  a/g  %  1  or  slightly  higher.  The  wavy  fin 
has  comer  angles  greater  than  90  degrees  whereas  the  corrugated  channels 
of  problem  CRACFLO  all  have  90  degree  bends.  Additionally,  the  comers 
of  the  wavy  fin  are  rounded  but  the  corrugated  channel  comers  are  sharp. 

These  differences  in  geometry  would  tend  to  reduce  the  size  of 
the  recirculation  zones  by  providing  a  more  streamlined  shape  with  less 
turning  of  the  flow.  The  rounding  of  the  comers  would  eliminate  the 
secondary  recirculation  zones  in  the  concave  comers  and  would  reduce 
the  separation  at  the  convex  comers.  Thus,  it  is  not  surprising  that 
Figure  4-3  does  not  resemble  the  upper  curve  of  Figure  4-12. 

Curve  Fits  and  Error  Analysis 

The  average  friction  factor  curves  have  been  correlated  for  the 
range  10  <_  Re  _<  1000.  The  resulting  curve  fit  is  a  third-order  poly¬ 
nomial  of  the  form: 

2  3 

f  -  cq  +  c^  log  Re  +  C2(log  Re)  +  c^Clog  Re) 


(4-5) 
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(4-9) 
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Equation  (4-9)  was  then  substituted  into  equation  (4-6)  to  estimate  the 
overall  accuracy.  Unfortunately,  the  accuracy  of  the  computed  results  is 
a  subjective  estimate  and  has  errors  of  its  own  which  presently  cannot  be 
determined. 

The  effect  of  a/g  of  f  Re  held  constant 

The  effect  of  varying  the  geometry  ratio,  while  holding  the  Rey¬ 
nolds  number  constant  at  Re  ■  250,  was  investigated  using  a  grid  having 
12  nodes  cross-stream.  This  spacing  places  the  numerical  value  of  the 
computer  results  in  question,  but  not  necessarily  the  relative  magnitudes. 

The  average  friction  factor  is  plotted  against  the  geometry  ratio 

in  Figure  4-13.  The  average  friction  factor  rises  to  a  maximum  at 

o. 

a/g  =  1  and  then  decreases  in  an  exponential  fashion  apparently  approach¬ 
ing  a  constant  value  as  a/g  increases.  It  was  previously  mentioned  that 
the  average  friction  factor  should  approach  the  friction  factor  for  laminar 
flow  between  parallel  plates  as  a  ■+■  ».  The  trend  toward  a  constant  value 
in  Figure  4-13  supports  this  hypothesis  but  it  appears  that  the  corner 
recirculation  continues  to  dominate  the  pressure  drop  even  at  a/g  £  10. 

General  Discussion 

The  results  presented  in  this  report  are  too  incomplete  to  be  of 
great  utility.  The  correlations  may  be  of  some  aid  but  only  for  the  two 
geometries  investigated.  However,  the  potential  of  program  CRACFLO  to 
serve  as  an  aid  in  the  engineering  design  of  corrugated  fin  heat  exchangers 
in  the  future  is  not,  by  any  means,  to  be  overlooked. 


Perhaps  Che  most  useful  observation  to  be  made  on  the  present 
results  is  that  laminar  flow  in  a  corrugated  channel  may  not  yield  the 
average  friction  factor  -  Reynolds  number  relationship 

f  *  Re-1  (4-10) 

This  observation  has  implications  when  one  attempts  to  predict  flow  pat¬ 
terns  in  a  large  rock  fracture  and  flow  through  a  crack  in  a  pressure 
vessel  [Button  et  al.  1978].  An  actual  crack  may  have  sharp  or  jagged 
walls  generating  recirculation  zones  which  affect  the  local  friction 
factor  as  in  a  corrugated  channel.  The  present  practice  of  assuming  a 
D'Arcy  relation  or  the  presence  of  smooth  walls  to  simplify  the  friction 
factor  calculation  is  expected  to  be  in  error  unless  the  Reynolds  num¬ 
ber  is  very  low. 

4.2  Problem  RUFF 

The  efforts  to  apply  the  IP  method  to  problem  RUFF  were  initi¬ 
ated  very  close  to  the  end  of  the  author’s  research.  The  direct  conse¬ 
quence  was  that  results  for  only  one  condition  -  geometry,  Prandtl  and 
Reynolds  numbers  and  thermal  conductivity  ratio  -  were  produced.  The 
discussion  here  will  be  brief  as  the  limited  results  do  not  warrant  a 
lengthy  discourse. 

Program  RUFF  determines  local  and  average  friction  factors  and 
Nusselt  numbers  in  addition  to  the  u  and  v  velocities,  pressures  and 
temperatures.  All  results  are  non-dimensional.  The  local  friction 
factors  and  Nusselt  numbers  are  determined  on  the  vertical  sides  of  the 
ribs  as  well  as  on  the  horizontal  surfaces. 
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Definitions 

The  local  friction  factor  is  defined  as 


f 


w 


(4-11) 


where  ft  is  a  unit  vector  perpendicular  to  the  surface,  u-  is  the  velocity 
vector,  and  the  subscript  "w"  indicates  the  values  are  determined  at  the 
fluid-solid  interface  or  wall. 

The  local  Nusselt  number  is  defined  as 


q  y 

kf  (Tw  -  V 


(4-12) 


where  is  the  local  heat  flux  at  the  wall,  Tw  is  the  local  wall  temper¬ 
ature,  is  the  greatest  distance  between  the  plate  wall  and  centerline 
(see  Figure  1-5) ,  k^  is  the  thermal  conductivity  of  the  fluid  and  the 
subscript  "B"  denotes  a  bulk  quantity. 

The  average  friction  factor  and  Nusselt  number  are  defined  such 
that  they  are  more  useful  for  general  design  considerations.  The  average 
friction  factor  is  defined  hy 


f 


2y_  A  P 

_  _£_  _ 

pUg7  (s  +  w) 


(4-13) 


where  AP  is  defined  in  equation  (4-2)  replacing  g  with  y„,  s  and  w 

cp  t 

are  the  rib  separation  and  rib  width  (see  Figure  1-5) ;  their  sum  is  a  rib 


cell  length. 
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The  average  Nusselt  number  is  defined  by 

q’cell  7a 


where 


Nu 


C 


kf(I»  -  tb  'l 

\  av  a vj 


=  jc  dl 

qcell  (s  +  w) 


fc  T 


di 


w  (s  +  w  +  2h) 
av 


(4-14) 


(4-15) 


(4-16) 


ts  +  w 


Tb  dx 


av 


(s  +  w) 


(4-17) 


The  integrals  for  ^  ^  and  T^  are  performed  along  the  contour  of  the 

av 

fluid  -  solid  interface  and  l  is  a  co-ordinate  for  distances  along  the 
contour. 

The  result  of  using  and  q'^^  in  these  average  quantities 

is  to  submerge  the  effect  of  the  rib  into  them.  This  allows  the  heat  ex¬ 
changer  designer  to  perform  his  job  without  concerning  himself  with 
details,  (which  are  available  should  they  be  of  interest).  One  should 
note  that  for  a  constant  heat  flux  applied  to  the  outside  of  the  repeated- 
rib  plate,  as  is  done  in  problem  RUFF,  that 


qcell  "  qP 


where  the  subscript  "p"  refers  to  the  plate. 


Additional  information  provided  by  program  RUFF  is  the  local  flow 


angles  and  the  maximum  temperature  in  a  plane  perpendicular  to  the  flow 
direction.  The  latter  allows  the  heat  exchanger  designer  to  locate  the 
hottest  points  in  the  heat  exchanger  and  compensate  as  required. 

The  Reynolds  number  used  in  program  RUFF  is  based  on  the  wall  to 
centerline  distance,  y  ,  because  this  is  the  characteristic  length  used 
to  non-dimensionalize  the  geometry 


Results 

A  single  solution  was  found  for  the  following  conditions 


Re  -  100 

V 

1.0 

Pr  =  0.7 

s  * 

0.25 

fluid:  air 

h  * 

0.04 

plate:  stainless  steel 

w  * 

0.12 

k.  ,k-  -  724 
p/  f 

0.21 

The  dimensions  are  non-dimensional  and  derived  from  the  actual  dimensions 
of  a  plate  with  repeated  ribs  to  be  used  in  an  experimental  test  rig  at 
Che  University  of  Arizona. 

A  heat  flux,  ,  applied  to  the  external  surface  of  the  plate  was 
modelled  as  well.  One  should  consult  chapter  two  for  a  discussion  on  the 
non-dimensionalization  of  the  energy  equation  and  an  explanation  of  how 
was  imposed. 

The  average  Nusselt  number  calculated  was  8.48  and  the  average 
friction  factor  calculated  was  0.07.  For  comparison,  the  Nusselt  number 


and  friction  factor  for  laminar  flow  between  parallel  plates  at  the  same 
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Reynolds  number  are  8.235  and  0.06,  respectively,  [Kays,  1966].  The  cal¬ 
culated  average  Nusselt  number  is  2.9%  greater  than  that  for  laminar  flow 
between  parallel  plates.  The  calculated  average  friction  factor  is  16.7% 
greater  than  the  friction  factor  for  laminar  flow  between  parallel  plates 
which  is  a  somewhat  unexpected  result.  According  to  Schlichting  [1968], 
the  roughness  element  investigated  should  have  little  effect  on  the  fric¬ 
tion  factor  in  laminar  flow.  This  might  be  attributable  to  a  too  permis¬ 
sive  overall  convergence  criterion  resulting  in  a  not  fully  converged 
solution.  The  convergence  criterion,  y  ,  was  0.005. 

The  local  Nusselt  numbers  are  plotted  in  Figure  4-14.  The 
abscissa  is  the  non-dimensional  distance  along  the  contour  of  the  fluid- 
plate  interface  measured  from  the  upstream  corresponding  plane.  One  may 
note  the  periodicity  of  the  results  and  that  the  local  Nusselt  number  is 
negative  at  two  points  indicating  heat  transfer  to  the  plate.  This  cor¬ 
responds  to  the  point  in  each  rib  indicated  by  the  circled  dot  in  Figure 
4-16. 

The  local  friction  factors  are  plotted  in  Figure  4-15.  The 
abscissa  is  the  same  as  was  used  in  Figure  4-14.  Here  also  one  can  note 
the  periodicity  of  the  results  and  the  negative  values  of  the  friction 
factor.  The  negative  values  occur  where  recirculation  is  present  and 
simply  implies  flow  in  the  negative  direction. 

The  maximum  temperature  at  all  cross-sections  occurs  at  the  ex¬ 
ternal  surface  of  the  plate  where  q"  is  applied.  This  result  is  reason- 

P 


able  and  predictable. 
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A  plot  of  the  local  flow  angles  is  presented  in  Figure  4-16  as  an 
aid  for  flow  visualization.  Only  the  region  above  the  rib  where  the  flow 
is  perturbed  noticeably  from  a  direction  parallel  with  the  centerline  is 
plotted.  A  small  recirculation  zone  occurs  before  and  after  the  rib, 
the  leeward  recirculation  zone  being  the  larger. 

The  velocity  profile  at  a  boundary  cross-section  is  plotted  in 
Figure  4-17.  The  co-ordinate  y^  is  the  distance  from  the  fluid-plate 
interface  measured  parallel  to  the  y-axis.  The  position  of  this  cross- 
section  relative  to  the  rib  can  be  reckoned  from  Figure  4-16  where  the 
left  or  right-hand  boundary  is  the  corresponding  plane  location.  A  lam¬ 
inar  parabolic  profile  is  plotted  as  a  solid  line  for  comparison.  The 
effect  of  the  rib  can  readily  be  seen  and  its  influence  extends  essential¬ 
ly  to  the  centerline.  Near  the  wall,  the  flow  is  still  recovering  from 
the  separation  caused  by  passing  over  the  rib. 
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CHAPTER  5 


CONCLUSIONS  AND  RECOMMENDATIONS 


5.1  Conclusions  for  Problem  CRACFLO 

A  computer  program  has  been  successfully  developed  for  numerical¬ 
ly  solving  the  fluid  flow  field  subject  to  spatially  periodic  boundary 
conditions  in  corrugated  channels  with  right  angle  bends.  Preliminary 
results  have  been  obtained. 

The  flow  in  this  corrugated  channel  is  seen  to  be  more  complex 
than  might  initially  be  anticipated;  it  exhibits  the  formation  of  a  large 
primary  recirculation  zone  and  a  small  secondary  recirculation  zone  as  the 
Reynolds  number  increases. 

As  a  consequence,  the  average  friction  factor  passes  through  a 
minimum  as  the  Reynolds  number  is  varied  while  the  geometry  ratio,  indi¬ 
cating  the  distance  between  successive  bends,  is  held  constant. 

Correlations  of  the  average  friction  factor  as  a  function  of  Rey¬ 
nolds  number  for  the  two  geometry  ratios  investigated,  0.0  and  0.5, 
were  determined.  The  form  of  the  fit  is 


2  3 

f  =  Cq  +  c^  log  Re  +  c2(log  Re)  +  c3(log  Re) 


where  the  values  of  the  constants  are 

a/g  =  0.0 

Cq  =  3.7107 

c1  =  -4.2743 
c2  =  1.7052 

c3  -  -0.22362 


a/g  3  0.5 

Cq  *  5.0625 

c  =  -6.633 
c,  =  3.2271 

c3  =  0.4665 
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The  accuracy  of  the  curve  fits  are  estimated  as  being  in  the  range  -4% 
to  -16%. 


Additionally,  the  effect  of  varying  the  geometry  ratio  while 
holding  the  Reynolds  number  constant  on  the  average  friction  factor  was 
investigated.  It  was  found  that  the  average  friction  factor  passed 
through  a  maximum  as  the  geometry  ratio  was  varied  and  approached  a  con¬ 
stant  value  as  the  geometry  ratio  grew  large.  The  latter  phenomena  was 
attributed  to  the  lessened  effect  of  the  comer  regions  containing  the 
recirculation  zones  on  the  entire  flow  in  duct  between  corresponding 
planes. 

Finally,  the  observed  behavior  of  the  average  friction  factor 
in  corrugated  channels  raises  doubts  on  the  applicability  of  using  the 
relation 

f  «  Re-1 

for  calculating  friction  factors  in  cracks  or  fractures  with  rough  walls. 

5.2  Recommendations  for  Problem  CRACFLO 

Recommendations  for  further  work  on  problem  CRACFLO  can  be 
grouped  into  three  categories:  a)  further  refinement  and  testing  of 
program  CRACFLO,  b)  production  of  more  extensive  predictions  for  laminar 
flow  in  corrugated  channels,  and  c)  application  of  program  CRACFLO  to  new 
problems . 

Program  CRACFLO  suffers  from  a  lack  of  versatility.  The  program 
should  be  modified  to  have  the  following  capabilities: 


J 
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1)  Ability  to  accurately  solve  for  flows  with  geometry 
ratios  to  25. 

2)  Ability  to  solve  the  thermal  energy  equation. 

3)  Ability  to  calculate  local  wall  shear  stresses  or  fric¬ 
tion  factors  and  local  Nusselt  numbers. 

4)  Ability  to  move  the  corresponding  planes  relative  to  the 
comers . 

5)  Ability  to  solve  for  transitional  and  fully  turbulent 
flows  as  well  as  laminar  flows. 

6)  Ability  to  perform  solutions  in  less  than  300  seconds 
on  a  CDC  CYBER  175. 

The  first  and  fourth  modifications  allow  for  further  testing  of 
program  CRACFLO,  particularly  to  investigate  the  behavior  of  the  average 
friction  factor  as  the  geometry  ratio  becomes  large.  The  remaining  modi¬ 
fications  are  to  increase  the  utility  of  program  CRACFLO  with  modifica¬ 
tion  six  perhaps  the  most  important  as  program  CRACFLO  presently  requires 
considerable  execution  time  and  is  therefore  expensive  to  use. 

An  important  test  of  program  CRACFLO  is  a  laboratory  experiment. 
This  experiment  is  crucial  for  establishing  the  validity  of  predictions 
made  by  program  CRACFLO.  Minimum  goals  of  the  experiment  would  be: 

a)  Determine  the  length  of  the  entrance  region  in  a  corru¬ 
gated  channel  prior  to  the  establishment  of  periodic 
flow, 

b)  Obtain  pressure  drop  data  to  allow  comparison  of  the 
average  friction  factors. 

c)  Perform  flow  visualization  to  chart  the  size  and  dev¬ 
elopment  of  the  recirculation  zones. 

The  limited  predictions  presented  in  this  report  should  be  ex¬ 
tended  and  new  predictions  produced.  Average  friction  factors  for  a 
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geometry  ratio  of  0.0  should  be  determined  for  1  <_  Re  <  10.  New  pre¬ 
dictions  of  the  average  friction  factor  over  the  range  1  <_  Re  <_  1000 
should  be  produced  for  geometry  ratios  of  1.0,  1.5  and  2.0.  Inspection 
of  the  geometries  investigated  by  Kays  and  London  [1964]  indicates  the 
geometry  ratio  range  0.0  _<_  a/g  <_  2.0  should  cover  most  of  the  geometries 
that  might  be  found  in  industry.  The  suggested  new  work  should  be  re¬ 
garded  as  a  minimum  goal  and  more  predictions  as  desirable. 

A  new  application  would  be  to  apply  the  results  of  program  CRACFLO 

as  input  to  a  program  predicting  laminar  flow  in  cracks.  This  application 
requires  a  correspondence  between  rock  fracture  geometries  and  the  geome¬ 
try  ratio  to  be  determined.  This  would  be  a  difficult  task,  but  has  poten¬ 
tially  great  rewards  for  predicting  the  performance  in  applications  such 
as  Hot-Dry-Rock  geothermal  reservoirs. 

A  second  new  application  is  to  modify  program  CRACFLO  to  predict 
the  performance  of  corrugated  absorber  plates  in  air-cooled  solar  collec¬ 
tors,  The  use  of  solar  collectors  for  domestic  heating  is  a  contemporary, 
but  not  a  mature  technology.  Hence,  there  exists  a  need  for  investigat¬ 
ing  new  methods  to  increase  collector  efficiency,  and/or  reduce  the  cost 
per  square  meter  or  unit  of  collected  energy.  A  corrugated  absorber  plate 
may  allow  one  or  both  of  these  goals  to  be  attained  by  improving  che  per¬ 
formance  by  heat  transfer  augmentation  [McEligot  and  Bankston,  1979]. 

5. 3  Conclusions  for  Problem  RUFF 

A  program  has  been  successfully  developed  for  numerically  solving 
for  laminar  flow  and  heat  transfer  subject  to  periodic  conditions  in 
parallel  plate  ducts  with  repeated  rib  roughness  and  a  constant  imposed 
heat  flux  applied  to  the  external  plate  surface. 


74 


A  single  case  was  investigated.  An  average  Nusselt  number  3% 
greater  and  an  average  friction  factor  16.7%  greater  than  those  for  lam¬ 
inar  flow  between  parallel  plates  were  found.  The  average  friction  fac¬ 
tor  result  is  believed  to  reflect  a  not  fully  converged  solution. 

Local  Nusselt  numbers  and  friction  factors  were  determined  as 

well.  The  results  indicate  that  heat  is  transferred  to  the  plate  at  the 
base  of  the  leeward  side  of  the  rib  and  that  recirculation  zones  occur 

before  and  after  the  rib. 

The  leeward  recirculation  zone  is  found  to  be  larger  than  the 
windward  recirculation  zone.  The  rib  affects  the  velocity  profile  at  a 
plane  approximately  midway  between  it  and  the  following  rib,  from  the 
wall  to  essentially  the  centerline. 

5 . 4  Recommendations  for  Problem  RUFF 

Recommendations  for  further  work  on  problem  RUFF  can  be  grouped 
into  two  categories:  a)  further  refinement  and  testing  of  program  RUFF, 
and  b)  application  of  program  RUFF  to  problems. 

Program  RUFF  was  brought  to  an  operational  state  by  the  author  in 
the  course  of  his  research,  but  many  tests  have  yet  to  be  performed.  The 
effect  on  the  calculated  flow  field  of  the  boundary  plane  location  rela¬ 
tive  to  the  rib  must  be  investigated.  Grid  independence  tests  must  be 
performed  and  the  effect  of  alternate  initial  flow  fields  investigated. 
The  overall  convergence  criterion  must  be  decreased  until  the  average 
friction  factor  becomes  insensitive  to  it. 

A  major  goal  is  to  reduce  the  computational  time  presently  re¬ 
quired  by  program  RUFF.  The  cyclic  tri-diagonal  matrix  method  mentioned 
in  Chapter  3  and  described  by  Faas  [1979;  Appendix  D]  may  provide  a  solu- 
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tion  to  this  problem.  Another  major  goal  is  to  modify  RUFF  to  handle 
transitional  and  turbulent  flows  as  these  are  of  greater  interest. 

The  existing  literature  should  be  surveyed  to  find  experimental 
data  to  aid  in  verifying  program  RUFF  and  aiding  in  its  development. 
Conducting  laboratory  experiments  may  be  necessary  as  well. 

Program  RUFF  should  eventually  be  used  in  the  design  and  opti¬ 
mization  of  heat  transfer  surface  using  repeated  rib  roughness  to  enhance 
heat  transfer.  The  current  grid  scheme  of  program  RUFF  restricts  appli¬ 
cation  so  only  small  roughness  elements  can  be  treated;  modifications 
should  be  made  to  allow  roughness  elements  on  the  order  of  half  the  plate 
spacing  to  be  investigated.  Other  geometries,  such  as  staggered  ribs  and 
ribs  on  a  single  wall  of  the  duct  are  extensions  of  program  RUFF  that 
should  be  easily  realized  and  of  interest.  Program  RUFF  may  have  appli¬ 
cations  in  predicting  flow  in  cracks  as  well  as  Program  CRACFLO. 
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Table  4-1.  Friction  factor  as  a  function 
of  Reynolds  number 
for  a/g  =  0.0  and  N^  =  24 

Non-Dimensional  ^cp 


Reynolds  Number, 

Re  Pressure, 

Friction  Factor 

10.0 

0.9181 

0.9181 

50.0 

0.2572 

0.2572 

100.0 

0.194 

0.194 

250.0 

0.183 

0.183 

500.0 

0.186 

0.186 

1000.0 

0.197 

0.197 

Table  4-2.  Friction  factor  as  a  function 
of  Reynolds  number 
for  a/g  *  0.5  and  Nq  *  24 

Reynolds  Number, 

Non-Dimensional 
Re  Pressure, 

IP 

CP 

OtJ,  *■ 

0 

Friction  Factor 

10,0 

1.786 

1.19 

50.0 

1.230 

0.82 

100.0 

1.540 

1.03 

250.0 

1,921 

1.28 

500.0 

2.181 

1.45 

1000.0 

2.41 

1.61 
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Table  4-3.  Friction  factor  as  a  function 
of  Reynolds  number 
for  a/g  »  0.5  and  N  =  12 


Reynolds  Number,  Re 

AP 

Non-Dimensional  co 

Pressure  Drop,  pU^ 

Friction  Fa< 

1.0 

15.37 

10.3 

2.0 

7.728 

5.15 

5.0 

3.203 

2.14 

10.0 

1.783 

1.19 

25.0 

1.135 

0.757 

50.0 

1.095 

0.730 

100.0 

1.255 

0.837 

250.0 

1.454 

0.969 

500.0 

1.656 

1.10 

1000.0 

1.861 

1.24 

78 


Table  4-4.  Friction  factor  as  a  function 
of  geometry  ratio, 
a/g,  for  Re  =*  250  and  *  12 


Geometry  Ratio,  a/g 

0.0 

0.5 

1.0 

1.5 

2.0 

5.0 

10.0 


Non-Dimensional  co 

Pressure  Drop,  pU, 2  Friction  Factor 

_ o  ______________ 


0.2857 

1.454 

2.165 

1.968 

1.608 

1.258 

1.601 


0.286 

0.969 

1.08 

0.787 

0.536 

0.210 

0.146 
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